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ABSTRACT 

Many  active  materials  exhibit  nonlinearities  and  hysteresis  when  driven  at  field  levels  necessary  to  meet 
stringent  performance  criteria  in  high  performance  applications.  This  often  requires  nonlinear  control  designs 
to  effectively  compensate  for  the  nonlinear,  hysteretic  field-coupled  material  behavior.  In  this  paper,  an  opti¬ 
mal  control  design  is  developed  to  accurately  track  a  reference  signal  using  magnetostrictive  transducers.  The 
methodology  can  be  directly  extended  to  transducers  employing  piezoelectric  materials  or  shape  memory  alloys 
(SMAs)  due  to  the  unified  nature  of  the  constitutive  model  employed  in  the  control  design.  The  constitutive 
model  is  based  on  a  framework  that  combines  energy  analysis  at  lattice  length  scales  with  stochastic  homoge¬ 
nizations  techniques  to  predict  macroscopic  material  behavior.  The  constitutive  model  is  incorporated  into  a 
finite  element  representation  of  the  magnetostrictive  transducer  which  provides  the  framework  for  developing 
the  finite-dimensional  nonlinear  control  design.  The  control  design  includes  an  open  loop  nonlinear  component 
computed  off-line  with  perturbation  feedback  around  the  optimal  state  trajectory.  Estimation  of  unmeasurable 
states  is  achieved  using  a  Kalman  filter.  The  hybrid  control  technique  provides  the  potential  for  real-time  control 
implementation  while  providing  robustness  with  regard  to  operating  uncertainties  and  unmodeled  dynamics. 

Keywords:  nonlinear  optimal  tracking,  magnetostrictive  actuators,  perturbation  control,  Kalman  filter 

1.  Introduction 

Smart  materials  and  structures  play  an  important  role  in  improving  performance  capabilities  of  aeronautic, 

aerospace,  automotive,  industrial,  biomedical,  and  military  systems.  The  need  for  high  power  density  devices 

that  can  operate  over  a  broad  frequency  range  with  precise  control  continues  to  draw  interests  in  transducer 
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designs  that  employ  smart  materials  such  as  piezoelectric  materials,  magnetostrictive  compounds  and  shape 
memory  alloys  (SMAs).  For  example,  piezoelectric  materials  offer  actuating  properties  with  nanoscale  reso¬ 
lution  and  large  actuation  forces  over  a  broad  range  of  frequencies  (Hz-MHz).  This  provides  advantages  in 
numerous  applications  such  as  high  resolution  acoustic  imaging  [1],  morphing  aircraft  control  surfaces  [2],  and 
high  resolution  nanopositioning  stages  [3,  4].  High  force  magnetostrictive  materials  such  as  Terfenol-D  provide 
broadband  capability  (DC  up  to  20  kHz)  with  forces  up  to  550  N  [5,  6],  which  are  utilized  in  many  industrial 
settings  such  as  high-speed  precision  machining  [7].  A  good  review  of  numerous  magnetostrictive  device  ap¬ 
plications  is  given  in  [8].  Shape  memory  alloys  operate  at  lower  frequencies  (<  100  Hz)  but  with  relatively 
large  power  density  which  provides  certain  advantages  in  applications  such  as  morphing  aircraft  structures  for 
improved  aerodynamic  performance  [9,  10]  and  reduction  of  jet  engine  noise  [11]. 

Whereas  smart  material  devices  offer  several  advantages  over  conventional  actuation  systems,  challenges 
associated  with  achieving  accurately  controlled  dynamics  require  careful  characterization  and  control  of  the 
field-coupled  material  behavior.  In  the  context  of  piezoelectric  and  magnetostrictive  materials,  a  small  input 
field  results  in  a  manifestation  of  local  polarization  or  magnetic  variants  moving  in  an  approximately  reversible 
manner  creating  macroscopic  strains,  which  are  approximately  linear.  Although  linear  models  and  control 
theory  often  provide  reasonable  approximations  when  the  input  fields  are  small,  moderate  to  large  fields  are 
often  necessary  to  meet  stringent  performance  criteria.  These  field  levels  typically  change  the  internal  state 
of  the  material  in  a  nonlinear  and  irreversible  manner  which  complicates  models  and  control  designs.  If  this 
behavior  is  neglected,  degradation  in  control  authority  may  result  from  unmodeled  constitutive  nonlinearities 
and  potential  instabilities  from  a  phase  lag  between  the  input  field  and  actuator  displacement.  It  is  often 
necessary  to  introduce  the  nonlinear  and  hysteretic  material  behavior  into  the  control  design  via  either  inverse 
filters  used  to  linearize  the  transducer  response  [7,  12]  or  direct  nonlinear  compensation  using  nonlinear  control 
designs  [13,  14,  15,  16].  Whereas  the  use  of  inverse  filters  can  improve  tracking  performance  by  effectively 
compensating  for  nonlinear  and  irreversible  material  behavior,  the  interpretation  of  the  transducer  input  is 
more  difficult  due  to  the  inversion  of  the  constitutive  law.  This  complexity  can  often  be  avoided  by  the  use  of 
nonlinear  control  design  which  circumvents  the  need  for  an  inverse  filter  and  directly  determines  the  nonlinear 
control  input.  The  inverse  filter  and  nonlinear  control  design  approaches  are  illustrated  in  Fig.  1. 

Model-based  control  designs  for  hysteretic  systems  are  not  new  and  have  received  considerable  attention 
in  areas  related  to  adaptive,  classical,  optimal  and  robust  control  [7,  12,  13,  14,  17,  18,  19].  These  models 
often  address  hysteretic  material  behavior  by  implementing  Preisach  operators  [12]  or  domain  wall  models 
[20].  Although  Preisach  models  can  reasonably  predict  nonlinear  and  hysteretic  material  behavior  observed 
in  smart  materials  such  as  piezoelectric,  magnetostrictive  and  shape  memory  alloys,  a  significant  number  of 


2 


Input  to  plant 


47 


(i) 


Hysteretic 

Transducer 


r 

7 

J 

) 

Hysteretic 

Transducer 

(ii) 


Figure  1:  (i)  Linear  control  design  employing  an  inverse  filter,  and  (ii)  nonlinear  control  design. 


non-physical  parameters  are  necessary  to  model  minor  loop  hysteresis.  Moreover,  significant  extensions  to  the 
classical  Preisach  theory  are  required  to  accommodate  relaxation,  reversible  behavior,  and  frequency-dependent 
dynamics.  Domain  wall  models  avoid  some  of  these  issues  by  incorporating  energetics  at  the  domain  length  scale 
to  reasonably  predict  unbiased  macroscopic  hysteretic  material  behavior.  However,  magnetostrictive  transducers 
are  typically  magnetized  by  a  permanent  magnet  to  provide  improved  actuation  properties  [5].  Similarly, 
piezoelectric  actuators  are  often  biased  with  a  static  electric  field  to  allow  increases  in  the  applied  unipolar 
field  and  corresponding  strain  without  inducing  reverse  ferroelectric  switching.  In  these  cases,  rate  independent 
biased  minor  loops  typically  exhibit  closure  properties  which  cannot  be  effectively  predicted  without  extensive 
modifications  of  the  domain  wall  model.  In  the  magnetization  model  developed  in  [21],  closure  of  biased  minor 
loops  was  enforced  by  a  priori  knowledge  of  the  turning  points  which  precludes  use  in  feedback  control  since 
the  states  are  not  known  in  advance.  When  inaccuracies  due  to  nonclosure  are  significant,  Preisach  models  or 
homogenized  energy  models  are  typically  employed  in  the  control  design. 

The  present  analysis  employs  the  homogenized  energy  framework  [22,  23,  24,  25]  in  the  control  design  to 
characterize  various  mechanisms  contributing  to  nonlinear  and  hysteretic  magnetostrictive  material  behavior. 
The  control  design  is  demonstrated  for  a  magnetostrictive  transducer  used  to  accurately  machine  an  out-of-round 
automotive  part  at  high  speed.  This  application  requires  a  100  pm  displacement  peak-to-peak  with  accuracy  of 
1  pm  -  2  pm  at  speeds  of  approximately  3000  rpm.  This  approach  is  similar  to  a  previous  investigation  of  high 
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speed  and  high  accuracy  machining  applications  [7]  that  utilized  an  inverse  filter  to  linearize  the  transducer 
response.  The  inverse  filter  was  implemented  in  a  robust  control  design  using  both  Hi  and  control  designs 
to  accurately  track  the  desired  profile.  The  control  design  developed  here  directly  determines  the  nonlinear 
control  input  by  implementing  the  magnetostrictive  constitutive  behavior  into  the  control  design.  Furthermore, 
the  control  input  is  determined  using  optimal  control  techniques  that  extend  recent  developments  in  hybrid 
nonlinear  optimal  vibration  control  of  plate  structures  [14].  The  previous  work  utilized  a  nonlinear  open  loop 
control  component  together  with  linear  perturbation  feedback  to  improve  robustness  to  operating  uncertainties 
in  actively  attenuating  plate  vibration  with  magnetostrictive  transducers.  The  control  design  presented  here 
required  modification  of  the  perturbation  feedback  control  when  tracking  a  reference  signal.  Additionally,  certain 
states  in  the  system  under  consideration  are  assumed  to  be  unmeasurable  which  motivated  implementing  a 
Kalman  filter  to  estimate  unobservable  states  in  the  perturbed  optimality  system.  This  approach  provides  a 
potential  method  for  real-time  state  feedback  of  perturbations  around  the  open  loop  nonlinear  control  for  high 
performance  tracking  applications. 

The  system  model  and  control  design  are  developed  as  follows.  First,  the  constitutive  model  is  briefly 
summarized.  A  structural  model  is  then  developed  using  one-dimensional  finite  elements  to  couple  the  magne¬ 
tostrictive  transducer  to  a  damped  oscillator  used  to  represent  boundary  conditions  between  the  transducer  and 
the  part  being  machined.  The  control  theory  is  then  developed.  Tracking  performance  is  demonstrated  using 
linear  optimal  control  theory  and  compared  to  the  nonlinear  optimal  control  approach.  Significant  improve¬ 
ments  in  accurate  position  control  are  obtained  when  the  magnetostrictive  constitutive  behavior  is  compensated 
using  the  open  loop  nonlinear  control.  To  improve  robustness  to  various  operating  uncertainties,  perturbation 
feedback  is  implemented  in  the  final  section.  A  linear  Kalman  filter  is  introduced  into  the  perturbed  optimality 
system  to  estimate  unobservable  states  in  the  transducer.  It  is  demonstrated  that  significant  improvements  in 
control  authority  are  achieved  when  using  the  perturbation  feedback  in  concert  with  the  Kalman  filter. 

2.  Magnetostrictive  Transducer  and  Hysteresis  Model 

Magnetostrictive  materials  deform  when  exposed  to  magnetic  fields  due  to  coupling  between  magnetic  mo¬ 
ment  orientations  and  elastic  interactions.  A  fundamental  challenge  in  effectively  utilizing  these  materials  in 
actuator  applications  deals  with  accurate  prediction  of  field-coupled  constitutive  nonlinearities  and  irreversibil¬ 
ities.  The  primary  mechanism  responsible  for  nonlinearities  and  hysteresis  is  domain  wall  motion  and  rotation 
around  pinning  sites  from  various  inhomogeneities  such  as  impurities,  inclusions  and  local  residual  stress  [26]. 
Whereas  this  behavior  can  be  mitigated  by  applying  only  small  fields,  moderate  to  high  input  fields  are  typi¬ 
cally  necessary  in  high  performance  applications.  This  requires  accommodating  constitutive  nonlinearities  and 
hysteresis  in  the  control  design  so  that  large  strain  and  high  force  capabilities  can  be  effectively  utilized. 
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Figure  2:  Terfenol-D  transducer  design  used  for  high  accuracy  machining  operations  of  cutting  out-of-round  objects. 


A  commonly  used  magnetostrictive  material  for  actuator  applications  is  Terfenol-D  (Tbo.3Dy0  7Fei. 9-1.95). 
This  composition  is  typically  employed  in  transducer  designs  due  to  its  broadband  capability  (DC  up  to  20  kHz), 
large  strain  (1,600  /re)  and  large  force  output  (550  N)  [27,  5].  This  material  is  often  manufactured  as  monolithic 
rods  for  incorporating  the  material  into  a  typical  transducer  design  as  depicted  in  Fig.  2.  The  transducer 
primarily  consists  of  the  Terfenol-D  rod,  a  surrounding  wound  wire  solenoid,  a  compression  bolt/spring  washer 
assembly,  and  a  permanent  magnet.  A  current  is  supplied  to  the  wound  wire  solenoid  which  produces  a  magnetic 
field  in  the  Terfenol-D  rod.  Magnetic  flux,  magnetization  and  strains  are  induced  by  the  magnetic  held.  The 
compression  bolt/spring  washer  assembly  prestresses  the  Terfenol-D  rod  to  avoid  tensile  loads  during  operation 
and  increased  actuator  displacement  [28].  The  permanent  magnet  serves  to  bias  the  magnetization  in  the 
Terfonol-D  rod  which  produces  uniform  flux  patterns  in  the  rod  and  bi-directional  strains  from  a  time- varying 
magnetic  held  input  with  zero  DC  offset.  This  avoids  the  ohmic  heating  associated  with  DC  bias  helds,  but 
increases  the  size  and  weight  of  the  transducer. 

2.1  Homogenized  Energy  Model 

The  homogenized  energy  model  is  based  on  a  multi-scale  modeling  approach  that  assumes  an  underlying 
distribution  of  material  property  relations  at  the  mesoscopic  length  scale  rather  than  assuming  uniform  ma¬ 
terial  properties  represented  by  material  constants.  A  stochastic  homogenization  technique  is  implemented  to 
determine  macroscopic  thermodynamic  parameters  from  the  distribution  of  the  local  material  properties.  The 
governing  equations  pertaining  to  the  control  design  are  summarized  here.  A  detailed  review  of  the  model  for 
ferromagnetic  materials  is  given  in  [22,  24]. 

We  focus  on  uniaxial  loading  of  rod-type  actuators  and  thus  assume  that  fields  and  material  coefficients  have 
been  reduced  to  scalar  coefficients  or  distributed  variables  in  the  direction  of  loading.  The  Gibbs  energy  at  the 
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mesoscopic  length  scale  is 


G{M ,  T)  =  (M,  T)  -  no HM  (1) 

where  'h  (M .  T)  is  the  Helmholtz  energy  detailed  in  [22],  no  is  the  permeability  of  free  space,  T  is  tempera¬ 
ture,  H  is  the  magnetic  field,  and  M  is  the  magnetization.  In  the  one-dimensional  case  considered  here,  the 
Helmholtz  energy  yields  a  double-well  potential  below  the  Curie  point  Tc  which  gives  rise  to  stable  spontaneous 
magnetization  with  equal  magnitude  in  the  positive  or  negative  directions. 

In  certain  cases,  thermal  activation  mechanisms  are  significant  and  must  be  included  in  the  constitutive 
model.  This  can  be  accomplished  by  introducing  the  Boltzman  relation 

h{G)  =  Ce~GV/kT  (2) 

which  quantifies  the  propability  /i  of  achieving  an  energy  level  G.  Here  kT /V  is  the  relative  thermal  energy 
where  V  is  a  representative  volume  element  at  the  mesoscopic  length  scale,  k  is  Bolztmann’s  constant,  and  the 
constant  C  is  specified  to  ensure  integration  to  unity.  From  a  physical  perspective,  large  thermal  energy  relative 
to  the  Gibbs  energy  reduces  the  barrier  for  magnetization  variants  to  switch  when  the  applied  field  is  near  the 
coercive  field.  This  reduces  the  sharp  transition  of  magnetic  variants  switching  from  the  positive  to  negative 
state  and  vice  versa. 

The  Boltzman  relation  is  utilized  to  determine  the  expected  magnetization  variants  governed  by 

/*  oo  /» —  Mi 

<Af+)  =  /  Mn(G)dM  ,  (M_)  =  /  Mn(G)dM  (3) 

J  Mi  J — oo 

where  ±Mj  are  the  positive  and  negative  inflection  points  in  the  Helmholtz  energy  definition. 

The  local  magnetization  is  subsequently  defined  by 

M  =  x+(M+)  +  x_{M_)  (4) 

where  x+  and  X-  respectively  denote  the  volume  fraction  of  magnetization  variants  having  positive  and  negative 
orientations.  The  differential  equations  governing  the  evolution  of  x+  and  X-  are  based  on  material-dependent 
parameters  that  define  the  likelihoods  that  moments  switch  from  positive  to  negative,  and  conversely.  These 
likelihoods  account  for  the  observed  thermal  relaxation  mechanisms.  Details  describing  the  governing  equations 
are  given  in  [22,  24]. 


6 


It  has  been  shown  [22]  that  in  cases  when  relative  thermal  effects  are  negligible  ( kT /V  is  small)  and  relaxation 
times  are  small  compared  to  the  drive  frequencies,  the  relation  Eq.  (4)  limits  to  the  piecewise  linear  relation 


M(H  +  Hf,Hc,Z)=Xm(H  +  HI)+6MR  (5) 

where  Xm  is  the  magnetic  susceptibility,  Mr  denotes  the  local  magnetization  variant  and  <5=1  for  local 
variants  aligned  in  the  positive  direction  and  8  =  —1  for  local  variants  aligned  in  the  negative  direction.  The 
interaction  field  Hj  represents  local  variations  in  the  field  due  to  material  inhomogeneities.  The  local  remanence 
magnetization  Mr  will  switch  when  the  diametrically  opposed  effective  field  ( He  =  H  +  Hj)  reaches  the  coercive 
field.  A  semi-infinite  set  of  magnetic  variants,  interaction  fields  and  coercive  fields  Hc  are  used  to  determine 
when  each  variant  switches.  Details  regarding  the  numerical  procedure  are  given  in  [22]. 

The  macroscopic  magnetization  is  computed  from  the  distribution  of  local  variants  via  the  relation 

nOO  nOO 

[ M(H )]  (t)  =  /  /  v(Hc,  Hj)  [M  {H  +  Hr,  Hc,  0]  {t)dHTdHc  (6) 

where  £  represents  the  initial  distribution  of  the  local  variants,  and  v{Hc,  Hi)  denotes  the  distribution  of  coercive 
and  interaction  fields.  To  facilitate  parameter  estimation,  the  following  form  of  the  distribution  is  used 


'(Hc,  Hi)  =  Cie-[ln(Hc/Hc)/2o]2e-^/262 


(7) 


where  Hc  is  the  average  coercive  field,  c  quantifies  the  coercive  field  variability,  b  is  the  variance  of  the  interaction 
field,  and  Ci  is  a  scaling  parameter.  Model  comparison  to  experimental  results  can  be  found  in  [22]. 

For  implementing  the  magnetostrictive  material  behavior  within  the  control  design,  the  forces  generated  by 
the  transducer  must  be  quantified.  This  is  provided  by  the  constitutive  law 


a  =  YMe-a1(M{H)-M0f  (8) 

representing  uniaxial  stress  in  the  Terfenol-D  rod  where  YM  is  the  elastic  modulus  at  constant  magnetization, 
£  is  the  linear  strain  component,  and  ai  is  the  magnetostrictive  coefficient.  It  is  assumed  that  stress  fields  are 
limited  to  the  linear  elastic  regime  where  ferroelastic  switching  is  negligible.  Here,  M{H)  is  the  magnetization 
computed  in  Eq.  (6)  whereas  Mo  includes  effects  from  the  permanent  magnet  and  the  initial  magnetized  state 
of  the  material. 

The  stress  computed  using  Eq.  (8)  includes  linear  stress-strain  behavior  as  well  as  nonlinear  and  hysteretic 
dependence  on  the  magnetic  field  through  the  M{H)  relation.  It  does  not  include  internal  damping  or  spatial 
dependence.  This  is  incorporated  in  the  following  section. 
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3.  Structural  Model 


To  facilitate  the  control  design,  the  constitutive  relations  given  by  Eqs.  (6)  and  (8)  are  used  to  develop 
a  system  model  that  quantifies  forces  and  displacements  when  a  magnetic  field  or  stress  is  applied  to  the 
magnetostrictive  transducer.  The  PDE  model  is  first  presented  and  then  formulated  as  an  ODE  system  through 
a  finite  element  discretization  in  space.  A  damped  oscillator  serves  as  boundary  conditions  at  the  end  of  the 
actuator  to  represent  the  structure  being  machined.  The  structural  model  is  illustrated  in  Fig.  3. 

A  balance  of  forces  for  the  structural  model  is  given  by  the  relation  [29] 


d2w  _  dNtot 
dt2  dx 


(9) 


where  the  density  of  the  actuator  is  denoted  by  p ,  the  cross-section  area  is  A  and  the  displacement  is  denoted 
by  w.  The  total  force  Ntot  acting  on  the  actuator  is  given  by  the  relation 


Ntot(t,  x)  =  YmA^  +  cDA^-t  +  Fmag{H)  +  Fd  (10) 

where  the  elastic  restoring  force  is  given  by  the  first  term  on  the  right  hand  side  of  the  equation  and  a  Kelvin- 
Voigt  damping  force  is  given  by  the  second  term.  The  linear  elastic  strain  component  is  defined  by  e  =  In 
comparison  to  Eq.  (8),  the  total  force  acting  on  the  Terfenol-D  rod  includes  the  additional  internal  damping 
force  and  the  external  disturbance  load  Fd.  The  coupling  force  Fmag  represents  forces  generated  by  an  applied 
magnetic  field  where 


Figure  3:  Magnetostrictive  transducer  with  damped  oscillator  used  to  quantify  loads  during  the  maching  operation. 
Disturbance  forces  along  the  actuator  are  given  by  Fd  and  the  control  input  is  u(t). 


As  illustrated  in  Fig.  3,  the  boundary  conditions  are  defined  by  a  zero  displacement  at  x  =  0  and  the  balance 


of  forces  at  x  =  L  given  by 


dw  d2w 

Ntot(t,  L)  =  - kLw(t,L )  -  cL  —  (t,  L)  -  mL-^p(t,L).  (12) 

The  initial  conditions  are  u>(0,x)  =  0  and  ^+(0,x)  =  0. 

The  strong  form  of  the  PDE  model  given  by  Eq.  (9)  can  be  written  in  the  variational  or  weak  form  for  finite 
element  implementation.  Multiplication  by  weight  functions  followed  by  integration  yields 


,  dw 


Y  A— — b  cdA 
ox 


d2w 

dxdt 


+  Fmag(H)  +  Fd 


kLw(t,  L)  +  cL  ^  (t,  L)  +  mL  (t,  L) 


(13) 


where  the  weight  functions  <fi(x)  are  required  to  satisfy  the  essential  boundary  condition  w(t,  0)  =  0  but  not  the 
natural  boundary  condition  given  by  Eq.  (12)  at  x  =  L.  The  weak  form  of  the  model  given  by  Eq.  (13)  is  used 
to  obtain  a  matrix  ODE  system. 


3.1  Approximation  Method 

An  approximate  solution  to  Eq.  (13)  is  obtained  by  discretizing  the  rod  in  space  followed  by  a  finite-difference 
approximation  of  the  resulting  temporal  system.  The  transducer  is  partitioned  into  equal  segments  over  the 
length  [0,  L\  with  points  xt  =  ih,  i  =  0, 1, . . . ,  N  and  step  size  h  =  L/N,  where  N  denotes  the  number  of  finite 
elements.  The  spatial  basis  {<pi}fL1  is  comprised  of  the  linear  interpolation  functions 


Mx)  =  *  < 

(X  -  Xi-i)  , 

Xi- 1  <  X  <  Xi 

(■ Xi+l  -  x)  , 

Xi  <  X  <  Xi+l 

(14) 

,  0 , 

otherwise 

for  i  =  1, . 

..,1V  —  1.  For  i  =  N,  the  basis  function  is  defined  to  be 

' — i  l-e 

II 

(£T 

% 

(x  -  XN-i)  , 

0  , 

Xat_i  <  X  <  Xn 

otherwise. 

(15) 

Details  describing  the  interpolation  functions  are  given  in  [22]. 

The  approximate  solution  to  Eq.  (13)  is  given  by  a  linear  superposition  of  the  interpolation  functions 


9 


(16) 


N 

wN{t,x)  =  5>(t)fc(*) 

1=1 

where  Wj(t)  are  the  nodal  displacement  solutions  along  the  length  of  the  transducer. 

An  ODE  matrix  system  is  obtained  by  substituting  wN(t,x)  into  Eq.  (13)  with  the  interpolation  functions 
employed  as  the  weight  functions.  This  gives  rise  to  a  second-order  matrix  equation 


Mw  +  Cw +  Kw  =  Fmag(H)b  +  fd  (17) 

where  the  displacement  solution  Wj(t)  has  been  written  in  vector  form  w(t)  =  [tui(f),  •  ■  • ,  ty/v(t)]  and  M  £  RNxN , 
C  £  RNxN  and  K  £  RNxN  denote  the  mass,  damping  and  stiffness  matrices.  The  vectors  b  £  RN  and  £  Rw 
include  the  integrated  interpolation  functions  related  to  the  control  input  and  disturbance  loads,  respectively. 
Components  contained  within  the  system  matrices  and  vectors  are  given  in  the  Appendix. 

Formulation  of  Eq.  (17)  as  a  first  order  system  yields 


x(t)  =  Ax(t)  +  [_B(u)](f)  +  G(t) 


z(0)  =  xq 


(18) 


where  x(t)  =  [w,w]T  should  not  to  be  confused  with  the  coordinate  x.  The  matrix  A  incorporates  the  mass, 
damping  and  stiffness  matrices  given  in  Eq.  (17)  and  [B(u)](t)  includes  the  nonlinear  input  where  u(t)  is  defined 
as  the  magnetic  field  to  be  consistent  with  the  control  literature.  External  disturbances  are  given  by  G.  The 
system  matrix  A,  input  vector  B{u )  and  disturbance  load  G{t)  are 


l 

O 

1=1 

0 

0 

A  = 

£ 

l 

o 

,  B(u)  =  Fmag(u) 

M-16 

,  G(t)  = 

M  -'Ut) 

The  identity  matrix,  with  dimension  N  x  N,  is  denoted  by  I. 

The  system  given  by  Eq.  (18)  must  be  discretized  in  time  for  simulation  purposes.  Here,  the  trapezoid  rule  is 
adopted  since  it  is  moderately  accurate,  A-stable,  and  requires  minimal  computer  storage  capacity.  A  standard 
trapezoidal  discretization  yields  the  iteration 


Xk+ i  —  Wit  +  [B{uk)  +  B(iik- t-i)  +  G(tk)  +  G(ffc+i)] 


Z(0)  =  Xq 


(20) 
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where  a  temporal  step  size  At  is  employed  giving  a  discretization  in  time  defined  by  tk  =  kAt.  The  values  Xk 
approximate  x(tk)-  The  matrices 


r  At  i 

-1 

r  At  j 

I - A 

I  H - A\ 

to 

to 

V  =  A  t 


(21) 


are  created  once  for  numerical  implementation  yielding  approximate  solutions  with  0(h 2,  (At)2)  accuracy. 


4.  Control  Design 

We  first  outline  the  general  tracking  problem  to  provide  relations  used  when  developing  the  nonlinear  optimal 
tracking  law.  The  tracking  problem  requires  defining  the  output  of  the  transducer  system.  This  is  developed 
by  assuming  that  the  position  of  the  cutting  tool  depicted  in  Fig.  (2)  can  be  measured  and  all  other  states  are 
unobservable.  The  output  to  Eq.  (18)  is  represented  by 


y(t)  =  Cx(t),  x  €  IR2JV 


(22) 


where  C  =  [0, . . . ,  0, 1, 0, . . . ,  0]  has  dimension  1  x  2 N.  The  cost  functional 


J  =  l(Cx{tf)  -  r(tf))TP(Cx(tf )  -  r{tf))  +  f  [H  -  A T{t)x(t))  dt  (23) 

Z  d  to 

is  defined  to  penalize  the  control  input  and  the  error  between  the  measurable  state  and  the  reference  signal 
r(t)  where  P  penalizes  large  terminal  values  on  the  observable  state,  H  is  the  Hamiltonian,  and  A (t)  is  a  set  of 
Lagrange  multipliers.  The  Hamiltonian  is, 


H  =  ^  \{Cx{t)  —  r(t))TQ(Cx(t)  —  r(f))  +  uT(t)Ru(t )] 

+AT  [Ax(t)  +  [B(u)](t)  +  G(t)} 


(24) 


where  penalities  on  the  observable  state  and  the  input  are  given  through  the  variables  Q  and  i?,  respectively. 

Through  the  construction  of  (7,  unobservable  states  are  not  included  in  the  performance  index.  The  cost 
functional  can  be  augmented  to  include  the  additional  states  without  conflict  in  performance  objectives  [30], 
but  this  is  not  considered  in  the  open  loop  control  design.  This  issue  is  treated  using  perturbation  feedback 
with  state  estimation  in  Section  4.3  to  minimize  potentially  undesirable  variations  in  unobservable  states. 
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To  determine  the  optimal  input,  the  minumum  of  Eq.  (23)  must  be  determined  under  the  constraint  of 
Eq.  (18)  which  is  incorporated  via  the  Lagrange  multipliers.  As  detailed  in  [31,  32],  the  stationary  condition  for 
the  Hamiltonian  yields  the  adjoint  relation 


A(i)  =  -AT\(t)  -  CTQCx{t)  +  CTQr(t) 


(25) 


and  control  condition 


u(t)  =  —R  1 


(26) 


The  boundary  conditions  are  applied  to  the  state  at  the  initial  time  and  the  costate  at  the  final  time 


x(t0)  =  x0 

A (tf)  =  CTP{Cx{tf)  -  r(tf)) . 


(27) 


These  equations  result  in  a  two-point  boundary  value  problem  that  is,  in  general,  challenging  to  solve  for  large 
systems.  In  the  following  section,  linear  control  theory  is  implemented  to  illustrate  that  poor  control  authority 
occurs  when  the  field  input  approaches  the  moderate  to  large  regime.  The  degradation  in  control  authority  is 
then  addressed  using  the  nonlinear  control  design.  It  is  demonstrated  that  by  directly  employing  the  constitutive 
model  in  the  optimality  system  and  solving  the  two-point  boundary  value  problem  through  a  combination  of 
analytical  and  numerical  techniques,  significant  improvements  in  tracking  a  reference  signal  are  obtained.  In  the 
last  section,  improvements  in  robustness  to  operating  uncertainties  are  realized  by  implementing  perturbation 
feedback  around  the  optimal  open  loop  trajectory  with  feedback  of  estimated  states  via  a  linear  Kalman  filter. 


4.1  Linear  Tracking  Control 

The  linear  control  design  requires  simplifying  the  nonlinear  input  given  by  the  coupling  force  in  Eq.  (11). 
The  force  is  linearized  by  taking  a  Taylor  series  expansion  about  the  fully  magnetized  state  of  the  material.  The 
control  input  in  this  case  is  limited  to  a  small  time  varying  input,  which  requires  only  a  consideration  of  first 
order  variations  in  the  magnetic  force. 

Linearization  about  Af0  yields  the  force  relation 


Fmag  —  2AuiAIqAM 


(28) 


where  AM  =  M  —  Mq. 
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A  linear  magnetic  constitutive  law  is  introduced  so  that  the  control  input  is  given  by  the  applied  magnetic 


field 


AM  =  M  -  M0  =  XmH 

where  H  is  the  change  in  magnetic  field  from  a  zero  DC  offset.  This  yields 


(29) 


Fmag  —  ~  d  Cl  1  -1  A)  A  m  ^ ' 

With  this  approximation,  the  corresponding  first-order  system  is 


(30) 


x(t)  =  Ax(t )  +  Bu(t)  +  G(t) 


x(0)  =  x0 


(31) 


where  u(t)  denotes  the  magnetic  field. 

The  state  constraint  in  Eq.  (31)  and  the  adjoint  condition  in  Eq.  (25)  yields  the  optimality  system 


x(t) 

_ 

A  ~BR~1Bt 

x(t) 

+ 

G(t ) 

m  _ 

-CTQC  -AT 

.  W  _ 

_  CTQr(t)  _ 

x(t0)  =  x0 


(32) 


x(tf)  =  Hfxitf)- 


Under  the  assumption  of  periodic,  steady-state  behavior,  we  consider  the  cost  functional 


J{u) 


'to 


( Cx(t )  —  r(t))  Q  ( Cx{t )  —  r(t))  +  uT(t)Ru(t) 


dt 


(33) 


where  r  is  the  largest  commensurate  frequency. 

Since  the  input  operator  B  is  linear,  a  fundamental  solution  matrix  is  determined  by  solving  the  algebraic 
Ricatti  equation 


AtU  +  UA  -  HBR-1BTn  +  CtQC  =  0 


(34) 


The  optimal  control  input  is  then  defined  by 


u*(t)  =  —R  1BT[Ux(t)  —  v(t)] 


(35) 
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where  the  variable  v(t)  £  R2Ar  is  the  solution  to  the  auxiliary  differential  equation 


v(t)  =  -  [A  -  BR~1BTn ] T  u(t)  -  CT Qr(t)  +  IIG(i) 
"(tf)  =  CTPr{tf) . 


(36) 


The  use  of  the  final  time  condition  results  from  approximation  of  the  periodic  boundary  conditions  to  facil¬ 
itate  implementation.  Details  regarding  this  approach  are  given  in  [33,  34].  The  control  trajectory  determined 
by  Eq.  (35)  is  used  in  the  following  numerical  examples  to  illustrate  the  need  for  nonlinear  control  when  the 
input  field  reaches  a  magnitude  that  induces  hysteresis. 


4.1.1  Numerical  Results:  Linear  Model 

We  demonstrate  tracking  control  using  the  linear  control  theory  to  illustrate  the  performance  in  tracking  a 
reference  signal  when  the  applied  field  level  is  either  small  or  moderate  to  large.  In  the  simulations,  we  apply  a 
large  input  field  to  the  rod  model  to  simulate  magnetization  of  the  Terfenol-D  rod  prior  to  control  simulation. 
The  set  of  parameters  given  in  Table  1  were  used  in  the  simulations.  The  values  for  the  penalities  on  the 
observable  state  and  input  were  Q  =  1  x  1012  and  R  =  1  x  10  10 .  The  step-size  used  in  all  simulations  was 
At  =  0.01  s. 

In  all  simulations,  the  disturbance  forces  Fd  were  set  to  zero.  It  was  determined  that  loads  in  excess  of  the 
coercive  stress  were  necessary  to  induce  noticeable  disturbance  vibrations  relative  to  the  tracking  displacement. 
This  eliminated  the  need  to  determine  the  optimal  input  for  disturbance  load  compensation  for  the  system 
under  consideration. 

Linear  tracking  control  is  demonstrated  using  a  reference  signal  that  starts  at  zero,  ramps  up,  oscillates  at 
a  specified  frequency  and  then  ramps  down  to  zero  in  accordance  with  the  machining  procedure  described  in 
the  Introduction.  The  control  is  turned  on  immediately  at  tQ  =  0.  In  Fig.  4,  the  reference  signal  is  limited  to 
small  displacements  so  that  the  constitutive  behavior  is  linear.  In  this  case,  excellent  tracking  is  achieved  as 
expected. 

When  the  reference  signal  is  increased  to  a  level  that  induces  hysteretic  magnetic  field-magnetization  be¬ 
havior,  degradation  in  control  authority  occurs.  Fig.  5(a)  illustrates  large  error  between  the  reference  signal 
and  the  controlled  displacement  which  clearly  does  not  achieve  the  1-2  /iin  tolerances  dictated  by  the  machin¬ 
ing  application  discussed  in  the  Introduction.  The  corresponding  hysteretic  H  —  M  behavior  is  illustrated  in 
Fig.  5(b). 
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(a)  (b) 

Figure  4:  (a)  Tracking  control  response  when  the  magnetostrictive  material  behavior  is  limited  to  the  linear  regime,  (b) 
The  H  —  M  behavior  was  computed  using  Eq.  (6). 


4.2  Nonlinear  Optimal  Tracking 

Accruate  tracking  in  the  presence  of  input  nonlinearities  and  hysteresis  is  addressed  in  this  section  by 
implementing  the  constitutive  behavior  directly  into  the  optimality  system  given  by  Eqs.  (18)  and  (25).  The 
methodology  follows  a  previous  approach  for  active  vibration  attenuation  of  beams  and  plates  [13,  14].  Key 

Table  1:  Parameters  used  in  the  tracking  control  problem. 


Magnetostrictive  and  Damped  Oscillator  Parameters 


kL  =  1.5  x  107  N/m 

Ym  =  3.0  x  1010  N/m 2 

cl  =  1.5  x  103  Ns/m 

A  =  5.0  x  10"4  in2 

n%L  =  0.47  kg 

L  =  0.1  m 

di  =  1.5  x  10-4  N/A2 

cr  =  3.7  x  106  Ns/m 

ci  =  6.1  x  10-5  m/A 

b  =  1.5  x  104  A/m 

Hc  =  3.3  x  103  A/m 

c=  0.65 

Mr  =  3.7  x  104  A/m 

Xm  =  15 

M0  =  2.75  x  105  A/m 
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Figure  5:  Linear  tracking  control  response  when  magnetostrictive  nonlinearities  and  hysteresis  are  present,  (a)  The 

controlled  response  is  given  by  (^— )  and  the  reference  signal  is  given  by  ( - ).  (b)  Corresponding  constitutive  behavior 

of  the  magnetostrictive  transducer  from  the  optimal  linear  control  input. 


equations  are  given  to  highlight  modifications  needed  to  track  a  reference  signal. 

In  the  nonlinear  model,  the  force  determined  from  Eq.  (28)  is  included  in  the  input  operator  [B(u)]  (t) 
where  the  force  generated  by  the  magnetic  field  has  been  linearized  about  the  biased  field  Mo-  Whereas  this 
approximation  neglects  the  quadratic  M  —  cr  relation  in  Eq.  (11),  it  includes  nonlinear  and  hysteretic  H  —  M 
relations.  Further  discussion  regarding  this  approximation  are  given  in  sections  4.2.1  and  4.3. 

The  resulting  optimality  system  is 

x{t)  Ax{t )  +  [B(u)](t)  +  G{t) 

A (t)  -AT\{t)  -  CTQCx(t)  +  CTQr(t) 

(37) 

x{t  o)  =  £o 

\{tf)  =  CTP{Cx{tf)-r{tf)). 

An  efficient  Ricatti  formulation  is  not  possible  in  this  case  due  to  the  nonlinear  nature  of  the  input  operator. 
The  two-point  boundary  value  problem  is  addressed  by  approximating  Eq.  (37)  or  the  equivalent  first-order 
system 
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z(t)  =  F(t,z ) 


E0z(t0)  =  [£o,0]t  (38) 

Efzitj)  =  [0,0]T 

where  z  =  [x(t),  X(t)]T . 

The  nonlinear  optimality  system  is  now  given  by 

Ax(t)  +  [B(u)](t)  +  G(t) 

F(t,z)  = 

-ATX(t)  -  CTQCx(t)  +  CTQr{t) 

(39) 

10  0  0 

Eq  —  ,  Ef  = 

0  0  -CTPC  I 

Here  I  denotes  an  identity  matrix  with  dimension  corresponding  to  the  number  of  basis  functions  employed  in 
the  spatial  approximation  of  the  state  variables. 

Various  methods  are  available  to  approximate  solutions  to  the  system  given  by  Eq.  (38)  such  as  finite 
differences  and  nonlinear  multiple  shooting  [35] .  The  finite  difference  approach  is  used  here  where  we  consider  a 
discretization  of  the  time  interval  [to,  tf]  with  a  uniform  mesh  having  stepsize  At  and  points  t0,  ti,  •  •  •  ,  tjv  =  t/- 
The  approximate  values  of  z  at  these  times  are  denoted  by  zq,  ■  ■  ■  ,  zjv-  A  central  difference  approximation  of 
the  temporal  derivative  then  yields  the  system 

^  izj+ 1  —  zj]  =  2  [-^(tj,  Zj)  +  F(tj_ |_i,  Zj+ 1)] 

E0z0  =  [y0,0]T  (40) 

EfZN  =  [0,  0] 

for  j  =  (],■■•  ,N—  1. 

Equation  (40)  can  be  expressed  as  the  problem  of  finding  Zh  =  [zo,  ■  •  •  ,  Zn ]  which  solves 

Hzh)  =  0  ■  (41) 

Equation  (41)  includes  the  optimality  system  at  each  time  step  and  the  boundary  conditions.  Details  are  given 
in  [13]. 

A  quasi-Newton  iteration  of  the  form 
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where  ££  solves 


~k+ 1  __  _fe  ,  fk 
Zji  I  S  h  5 


(42) 

(43) 


is  then  used  to  approximate  the  solution  to  the  nonlinear  system  given  by  Eq.  (41).  The  Jacobian  ^F'(z^)  has 
the  form 


where 


F'(zh) 


Sn-i 


(44) 


0 

I 


^4  fxB[u*] 

-CTQC  ~AT 


(45) 


for  Si .  The  representation  for  Ri  is  similar. 

Inversion  of  the  Jacobian  is  required  to  obtain  updates  to  the  states  through  the  iteration  procedure.  Direct 
inversion  is  not  feasible  when  a  large  number  of  basis  functions  are  required  to  discretize  the  structural  model. 
In  the  one-dimensional  rod  model,  this  is  typically  not  an  issue,  although  to  reduce  potential  numerical  errors 
in  inverting  the  Jacobian  an  analytic  LU  decomposition  is  employed.  The  form  of  the  equation  is  identical  to 
that  used  in  previous  investigations  [13,  14]  and  is  given  here  in  the  Appendix  for  convenience.  This  gives  rise 
to  the  following  representation  of  the  Jacobian 


T\zkh)  =  LU  (46) 

where  L  and  U  respectively  denote  lower  and  upper  triangular  matrices. 

The  solution  of  the  system  defined  by  Eq.  (43)  is  then  obtained  through  direct  solution  of  the  lower  triangular 
system  L( =  —R(zfc)  followed  by  direct  solution  of  the  upper  triangular  system  U =  (£.  This  leads  to  updates 
to  the  state  and  adjoint  soutions  z%+1. 
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(a)  (b) 


Figure  6:  Tracking  control  response  using  the  nonlinear  control  law.  (a)  The  controlled  response  (in  gm)  (^— )  and 

the  error  between  the  reference  signal  and  controlled  response  (in  fim)  ( - ).  (b)  Corresponding  H  —  M  constitutive 

behavior. 


4.2.1  Nonlinear  Simulation  Results 

The  nonlinear  optimal  tracking  design  provides  enhanced  control  authority  when  nonlinearities  and  hysteresis 
are  present.  This  is  demonstrated  by  prescribing  the  same  reference  signal  that  induced  hysteresis  in  Fig.  5(b) 
and  poor  tracking  in  Fig.  5(a).  The  transducer  nonlinearities  are  effectively  accommodated  using  the  nonlinear 
optimal  control  method  as  shown  in  Fig.  6.  In  this  simulation,  the  optimal  control  input  was  computed  assuming 
that  the  actuation  force  was  proportional  to  the  magnetization  as  previously  given  by  Eq.  (28).  In  Fig.  6,  the 
magnetostrictive  actuator  displacement  was  determined  from  the  control  input  using  the  quadratic  magnetization 
coupling  as  given  by  Eq.  (11).  This  illustrates  that  the  error  from  unmodeled  quadratic  coupling  is  minimal  for 
a  100  /mi  reference  trajectory.  The  same  penalties  on  the  observable  state  and  control  input  from  the  linear 
simulations  have  been  used  with  the  additional  penalty  on  the  final  state  given  by  P  =  CHCT . 

Remark  1:  In  the  present  model,  nonlinearities  were  included  in  the  input  operator  contained  within  the 
optimality  system,  Eq.  (37).  It  should  be  noted  that  this  method  neglects  nonlinearities  contained  within  the 
components  of  the  Jacobian  (Eq.  (45))  and  the  optimal  input  equation  (Eq.  (26)).  The  analytic  form  of  the  term 
-§^B  [it*]  in  Eq.  (45)  requires  explicit  knowledge  fo  the  Lagrange  multipliers  and  the  input  which  is  unknown 
in  the  nonlinear  case.  This  component  is  assumed  to  be  constant  which  results  in  a  suboptimal  controller  with 
Si  =  S  and  Ri  =  R  for  each  time  step.  It  should  also  be  noted  that  optimal  input  contains  nonlinearities  in 
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the  term  defined  by  dBg^u'>  ■  The  nonlinear  term  requires  computing  the  tangential  magnetic  susceptibility, 
Since  the  input  only  induces  minor  hysteresis  in  the  present  analysis,  changes  in  the  magnetic  susceptibility 
are  expected  to  be  small  and  are  neglected  in  the  nonlinear  control  simulations.  The  tangential  magnetic 
susceptiblity  is  assumed  to  be  equivalent  to  the  linear  coefficient  Xm  = 

Remark  2:  The  control  design  presented  here  has  focused  on  open  loop  nonlinear  control  that  is  computed 
off-line.  It  is  well  known  that  uncertainties  in  operating  conditions  or  unmodeled  dynamics  can  have  deleterious 
effects  on  control  performance.  This  is  addressed  in  the  following  section  using  perturbation  feedback.  The 
approach  is  similar  to  a  hybrid  nonlinear  optimal  control  design  for  attenuating  plate  vibration  [14].  While 
there  are  similarities  to  the  tracking  problem,  fundamental  differences  in  the  auxiliary  equation  (Eq.  (36))  are 
addressed  to  accurately  track  the  prescribed  reference  signal  when  operating  uncertainties  are  present. 

4.3  Perturbation  Control  and  Estimation 

In  the  previous  section,  nonlinear  open  loop  control  was  shown  to  provide  excellent  tracking  behavior  by 
compensating  for  nonlinear  and  hysteretic  field-coupled  magnetostrictive  material  behavior.  However,  any 
unmodeled  dynamics  not  accounted  for  in  the  constitutive  model  or  damped  oscillator  are  a  source  of  operating 
uncertainty  which  can  be  detrimental  in  obtaining  high  accuracy  tracking.  Moreover,  external  disturbances  or 
slight  delays  in  the  control  electronics  can  degrade  control  authority.  These  issues  are  addressed  by  implementing 
perturbation  feedback  around  the  optimal  open  loop  trajectory  computed  in  the  preceding  section. 

As  previously  noted,  the  system  under  consideration  contains  unmeasurable  states  defined  by  the  output 
equation  (Eq.  (22))  since  it  has  been  assumed  that  only  the  displacement  at  the  end  of  the  cutting  tool  is  known. 
State  feedback  requires  availability  of  all  the  states  which  necessitates  the  estimate  of  the  unmeasurable  states. 
This  is  accomplished  by  employing  a  Kalman  filter.  The  estimated  perturbations  are  then  used  to  determine 
first-order  variations  in  the  optimal  control  input.  The  development  of  the  perturbation  controller  and  estimator 
forms  the  well-known  linear  quadratic  Gaussian  (LQG)  control  problem  where  incomplete  state  information  is 
used  to  determine  the  LQ  regulator  described  in  Section  4.1  except  the  control  is  determined  from  feedback 
from  the  state  estimates  using  the  Kalman  filter.  Key  equations  used  in  developing  the  perturbation  controller 
and  Kalman  filter  are  given  here.  Additional  details  can  be  found  in  [14,  36,  37]. 

The  perturbed  system  consists  of  first-order  variations  in  the  constraint  Eq.  (18),  output  Eq.  (22)  and  initial 
conditions, 
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Sx(t)  =  ASx(t)  +  BSu(t)  +  SG(t) 


Sy(t)  =  C5x(t)  +  w(t)  (47) 

5x(t0)  =  Sx  o, 

where  Su,  Sx  and  Sy  are  first-order  variations  about  it*,  x*  and  y* .  The  nonlinear  input  operator  B(u)  has 
been  linearized  around  the  optimal  input  u*(t).  The  perturbation  in  plant  disturbances  is  denoted  by  SG(t). 
In  comparison  to  the  output  in  Eq.  (22),  measurement  noise  w(t)  has  been  included  in  the  first-order  variation 
of  the  output  to  illustrate  robustness  in  the  resulting  design.  Both  the  variation  in  system  noise  SG(t)  and  the 
measurement  noise  w(t)  are  assumed  to  be  random  processes  with  zero  mean.  The  covariance  for  the  system 
and  measurement  noise  is  defined  by  E  [<5G(t)<5GT(r)]  =  VS(t  —  r)  and  E  [w(t)wT (r)]  =  WS(t  —  r)  with  the 
joint  property  E  [<5G(f)ieT(r)]  =  0. 

The  cost  functional  given  by  Eq.  (23)  is  used  to  determine  the  optimal  perturbation  feedback.  Since  the  opti¬ 
mal  control  pair  x*(t))  minimizes  SJ,  the  second  order  terms  are  expanded  to  determine  (Su*(t),Sx*(t)). 

The  second  order  variation  in  J  is  given  by 


S2J  = 


f  [SxT  SuT ] 

Q  0 

Sx 

'to 

0  R 

Su 

dt 


(48) 


where  penalties  on  the  perturbed  states  and  input  are  given  by  Q  and  R,  respectively.  Here,  Q  £  jpj2Wx2W 
is  a  matrix  that  weights  both  the  displacement  and  velocity  at  the  end  of  the  cutting  tool.  This  minimizes 
the  potential  for  large  variations  in  the  cutting  tool  velocity.  A  similar  approach  is  given  in  [30].  The  form  of 
SxTQ5x  is  given  by 


SxTQSx  =  q\Sx2N  +  Q2Sx\n  (49) 

where  Sxn  is  the  variation  in  displacement  at  the  tip  of  the  cutting  tool  and  Sx^n  is  the  variation  in  velocity 
at  the  tip  of  the  cutting  tool.  The  parameters  q\  and  q-2  are  the  respective  penalties. 

Due  to  the  form  of  the  performance  index  given  by  Eq.  (48) ,  previously  discussed  LQR  theory  can  be  directly 
employed  to  obtain  6u*(t )  and  Sx*(t).  The  overall  control  for  the  system  is  then  taken  to  be  u*(t)  +Su*(t)  with 
the  optimal  state  given  by  x*(t)  +  Sx*(t).  The  perturbed  feedback  signal  is  given  by 

Su*(t)  =  —  R~1Bt  [n&r(£)  —  Su(t)]  (50) 
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where  II  is  the  solution  of  the  algebraic  Ricatti  equation  for  the  perturbed  system  previously  given  by  Eq.  (34). 
The  variable  8x(t)  is  the  estimate  of  the  first-order  variations  in  the  states  given  by  Eq.  (47)  which  will  be 
subsequently  determined  from  the  Kalman  filter. 

The  variable  Sv{i)  given  in  Eq.  (50)  is  a  feedforward  term  used  to  compensate  for  phase  delays  or  periodic 
external  disturbances.  It  is  governed  by  the  auxilliary  equation 


8v{t)  =  -  [A  -  BR~1BtU]T  Su(t)  -  CTQSr(t )  +  USG{t) 
<M tf)  =  CTPSr(tf) . 


(51) 


where  the  final  time  conditions  are  zero  in  the  present  machining  example. 

To  solve  Eq.  (51),  the  perturbation  in  the  disturbance  SG(t)  and  reference  signal  Sr(t )  must  be  known  or 
estimated.  The  perturbation  in  the  reference  signal  is  defined  by  Sr{t)  =  y(t)  —  r[t  +  St)  where  the  time  shift 
is  denoted  by  St.  As  mentioned  previously,  disturbance  loads  have  a  negligible  effect  on  tracking  performance, 
therefore  perturbations  in  the  disturbance  loads  SG{t)  are  not  considered  in  the  simulations. 

The  perturbation  feedback  control  given  by  Eq.  (50)  requires  estimation  of  the  states  in  Eq.  (47)  which  are 
given  by  the  relation 


Sx{t)  =  ASx(t)  +  BSuit)  +  L(t)  [6y(t)  —  CSx(t)] 
Sx(to)  =  Sxg 


(52) 


where  Sx o  is  the  estimated  initial  conditions  and  L(t)  is  the  Kalman  gain. 

The  Kalman  gain  is  determined  from  the  dynamics  of  the  error  between  the  first-order  variation  in  the  actual 
states  and  the  estimated  states.  Key  equations  used  in  determining  L(t)  are  given  here.  Additional  details  can 
be  found  in  [36,  37]. 

The  Kalman  gain  is  determined  from  the  error  covariance  T (t)  =  E  [e(t)eT(r)]  where  the  error  is  determined 
from  the  dynamic  relation 


e(t)  =  5x(t)  —  Sx(t).  (53) 

By  substituting  Eqs.  (47)  and  (52)  into  Eq.  (53)  the  error  is 

e(£)  =  $(£,  t0)e(t0)  +  f  4>(t,  r)  [<5G(r)  -  L(t)w(t)\  dr  (54) 

Jt0 

where  <!>(£,  to)  is  the  state  transition  matrix.  Through  manipulation  of  T (t)  outlined  in  [36,  37]  the  rate  of 
change  of  the  error  covariance  leads  to  the  differential  matrix  Ricatti  equation 
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f  (t)  =  AT(t)  +  T (t)AT  -  T(t)C'TW/"1C'T(i)  +  V 


(55) 


with  the  initial  conditions  Y(to)  =  To-  The  Kalman  gain  is  found  by  maximizing  the  negative  rate  of  change  of 
T  (t)  which  leads  to 


L(t)  =  T(t)(7TW_1.  (56) 

In  the  simulations  presented  here,  a  steady-state  Kalman  gain  L(t)  =  L  is  determined  by  solving  the  algebraic 
form  of  Eq.  (55)  where  T (t)  =  0.  This  leads  to  an  approximately  optimal  system  that  is  easier  to  implement 
numerically. 

In  summary,  the  matrix  system  governing  the  dynamics  of  the  perturbed  optimality  system 


Sx(t ) 

_ 

A  -BR-1BTn 

Sx(t) 

+ 

BR~1Bt 

8v(t)  + 

SG(t) 

Sx(t) 

LC  A  -  BR~lBTn-  LC 

§x(t) 

br_1bt 

Lw(t) 

(57) 

Sx(to)  =  Sx  o 
Sx(t-o)  =  Sxq. 

is  solved  to  obtain  the  system  response  when  perturbation  feedback  of  the  estimated  states  is  implemented. 

4.3.1  Numerical  Results 

The  effectiveness  of  the  perturbation  controller  and  estimator  is  illustrated  by  applying  an  initial  phase 
delay  of  St  =  0.03  s  to  the  Terfenol-D  transducer  to  model  the  effects  of  delays  in  the  electronics.  The  tracking 
performance  is  compared  to  the  cases  when  measurement  noise  is  either  negligible  or  significant.  In  these 
simulations,  the  system  response  was  assessed  when  the  measurement  noise  was  corrupted  by  a  signal  with  a 
100  //in  amplitude  and  compared  to  a  baseline  case  when  w(t)  =  0.  As  previously  noted,  whereas  disturbance 
loads  are  included  in  the  governing  equations,  these  forces  were  neglected  in  the  simulations  since  large  stresses 
in  excess  of  the  coercive  stress  were  necessary  to  induce  undesirable  vibration  in  the  rod.  A  similar  approach  is 
taken  here  where  perturbations  in  the  disturbance  forces  are  neglected.  The  quantities  used  for  the  penalties  in 
the  second  variation  of  the  cost  functional  were  q\  =  1  x  1016,  <72  =  1  x  1011  and  R  =  lx  10-12.  The  Kalman 
filter  gain  was  computed  using  a  value  of  V  =  1  x  108  as  a  design  parameter  for  the  covariance  of  6G(t)  and 
W  =  lx  10-7  for  the  covariance  of  w(t). 
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The  system  performance  with  and  without  perturbation  feedback  is  illustrated  in  Figs.  7  and  8,  respectively. 
As  expected,  the  open  loop  tracking  control  in  Fig.  7  results  in  an  offset  in  displacement  shifted  by  0.03  s  relative 
to  the  desired  reference  trajectory  leading  to  large  errors  in  tracking  the  reference  signal.  Significant  improve¬ 
ments  are  made  when  perturbation  feedback  control  is  employed.  Figure  8  compares  tracking  performance  when 
the  measurement  noise  is  zero  (Fig.  8(a))  to  performance  when  the  noise  is  nonzero  (Fig.  8(b)).  Whereas  the 
error  is  initially  large  from  the  phase  shift,  it  subsequently  reduces  to  values  comparable  to  the  open  loop  case 
when  delays  in  the  electronics  are  absence  (Fig.  6)  except  at  points  where  abrupt  transitions  in  displacement 
occur.  When  measurement  noise  is  present,  the  Kalman  filter  limits  the  error  to  approximately  the  amplitude 
of  w(t)  as  illustrated  in  Fig.  8(b). 


Figure  7:  The  effect  of  an  initial  time  delay  in  the  open  loop  control  signal  is  illustrated  where  (^— )  is  the  controlled 
response  and  ( - )  is  the  reference  signal. 
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10 


(a)  ( b ) 


Figure  8:  Perturbation  feedback  using  the  estimated  states  from  the  Kalman  filter:  (a)  measurement  noise  is  zero,  and 
(b)  measurement  noise  is  a  random  signal  with  an  amplitude  of  1  pm. 


5.  Concluding  Remarks 

A  model-based  nonlinear  optimal  control  design  was  developed  to  accurately  track  a  reference  signal  when 
nonlinear,  hysteretic  magnetostrictive  material  behavior  is  present.  Whereas  linear  control  was  effective  for 
tracking  the  desired  reference  signal  for  small  displacements,  significant  errors  occur  at  moderate  to  high  field 
inputs  if  the  constitutive  behavior  is  not  properly  accommodated  .  This  was  addressed  by  incorporating  the 
homogenized  energy  model  into  the  control  design  to  compensate  for  nonlinear,  hysteretic  field-coupled  mag¬ 
netostrictive  constitutive  behavior.  The  open  loop  nonlinear  control  method  provided  excellent  tracking  per¬ 
formance  at  field  levels  where  nonlinearities  and  hysteresis  are  present.  Although  this  enhanced  performance, 
open  loop  control  is  not  practical  in  applications  due  to  uncertainties  in  operating  conditions.  This  issue  was 
addressed  by  complementing  the  open  loop  nonlinear  control  with  linear  perturbation  feedback  using  LQG 
control.  This  hybrid  control  technique  included  a  Kalman  filter  in  the  perturbed  optimality  system  to  feedback 
estimations  of  variations  in  states  around  the  optimal  open  loop  trajectory.  This  approach  provides  a  potential 
method  for  real-time  control  when  certain  states  are  unobservable  by  computing  and  storing  the  open  loop 
nonlinear  control  and  applying  real-time  linear  feedback  of  the  perturbed  states. 

The  nonlinear  tracking  control  design  has  focused  on  compensating  for  rate-independent,  nonlinear  and 
hysteretic  magnetostrictive  material  behavior,  although  the  control  design  is  general  for  a  broad  class  of  ferroic 
materials  due  to  the  unified  properties  on  the  homogenized  energy  model.  The  constitutive  model  has  been 
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successful  in  modeling  the  behavior  of  ferroelectric  materials  and  shape  memory  alloys  [25].  Moreover,  the 
homogenized  energy  framework  is  applicable  to  rate-dependent  hysteresis  in  materials  where  thermal  relaxation 
mechanisms  are  present  or  the  drive  frequency  affects  the  level  of  hysteretic  losses.  Implementing  this  behavior 
into  nonlinear  control  designs  is  under  current  investigation. 
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Appendix 


A.l  Matrix  and  Vector  Relations 

The  linear  interpolation  functions  used  to  formulate  the  matrices  and  vectors  given  in  Section  3.1  are  given 
here.  By  substituting  the  approximate  displacements  in  Eq.  (16)  into  the  variational  form  of  the  dynamic 
equation,  Eq.  (13)  has  the  form 

JV  N  N  „L  [.L 

y^u)j(f)  /  pA(j>i(j)jdx  +  Wj ( t )  /  cpAfi^dx  +  ^  Wj(t)  /  YM Acjj'^dx  =  Fmag(H)  /  <$>\dx 

j=l  J0  j= 1  Jo  j=1  Jo  Jo 

+  /  Fd4>idx  -  ( kLwN(t)(j)N(L )  +  cLwN(t)4>N(L)  +  mLwN(j)N(L ))  4>N(L) . 

Jo 

By  comparing  the  above  relation  to  Eq.  (17),  the  global  mass,  stiffness,  and  damping  matrices  contain  the 
components 


[M] 


ij 


pA(t>i(j)jdx  , 


/  pA4>i(j)jdx  +  niL  , 

Jo 


i  N  or  j  N 
i  =  N  or  j  =  N 


J  cpA^jdx  , 


cpA^cji'jdx  +  cl 


i  N  or  j  N 


i  =  N  or  j  =  N 
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Ym  Afttfjdx  ,  i^N  or  j^N 

YM A^cpjdx  +  Ul  ,  *  =  N  or  j  =  N 

whereas  the  force  vectors  have  the  components 


r  L  _  f'L/ 

fd=  Fd(f>idx  ,  b  =  /  (p'j/lx. 

L  Jo  1  J*  Jo 

A. 2  LU  Decomposition 

The  Jacobian  given  by  Eq.  (44)  is  written  in  terms  of  a  LU  decomposition  to  avoid  matrix  inversion  when 
considering  large  systems.  The  analytic  form  of  the  decomposition  used  in  the  simulations  is 


So'Ro 

I 


I 


This  form  of  the  matrix  system  reduces  the  computation  to  solving  2N  x  2 N  systems  for  each  time  step 
separately.  This  provides  the  ability  to  solve  matrix  systems  containing  a  large  number  of  degrees  of  freedom 
while  maintaining  a  small  temporal  step  size. 
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